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1 | INTRODUCTION 


Protected areas are essential to help safeguard biodiversity from an- 
thropogenic pressures (Watson, Dudley, Segan, & Hockings, 2014). 
Since resources are limited, systematic analyses are used to reveal 
shortfalls in existing protected area systems (termed gap analyses) 
and develop cost-effective plans for expanding them (termed prior- 
itizations,; Margules & Pressey, 2000). Protected area systems can 
promote the long-term persistence of species by securing habitats 
and populations that are important for maintaining their evolu- 
tionary processes (Moritz, 2002). Specifically, neutral evolutionary 
processes create and regulate genetic diversity, and adaptive evo- 
lutionary processes enable species to overcome changes to their 
surrounding environments (reviewed in Funk, McKay, Hohenlohe, 
& Allendorf, 2012). As a consequence, there has been increasing 
interest in incorporating evolutionary processes into conservation 
planning (Carvalho, Brito, Crespo, & Possingham, 2011; Carvalho 
et al., 2017; Hanson, Rhodes, Possingham, & Fuller, 2018; Nielsen, 
Beger, Henriques, Selkoe, & von der Heyden, 2017; Paz-Vinas 
et al., 2018; Thomassen et al., 2011). 

Protected areas should ideally maintain and enhance adaptive 
and neutral evolutionary processes (Sgro, Lowe, & Hoffmann, 2011). 
Previous prioritizations have incorporated such processes at the in- 
traspecific level using maps of related attributes (hereafter, evolu- 
tionary attributes). For example, maps of neutral genetic diversity 
and gene flow have served as evolutionary attributes for neutral 
processes (Hanson, Fuller, & Rhodes, 2019; Hanson et al., 2018; 
Nielsen et al., 2017). On the other hand, maps of climatic refugia, 
environmental conditions, morphological variation and putatively 
adaptive genetic diversity have served as evolutionary attributes 
for adaptive processes (Cowling, Pressey, Rouget, & Lombard, 2003; 
Game, Lipsett-Moore, Saxon, Peterson, & Sheppard, 2011; Hanson, 
Rhodes, Riginos, & Fuller, 2017; Thomassen et al., 2011). Although 
an optimally sited prioritization would represent a comprehensive 
range of these evolutionary attributes for each species of interest 
(Beger et al., 2014), long-standing challenges in mapping and opera- 
tionalizing evolutionary processes have severely limited the scope of 
previous prioritizations. 

Here, we introduce a new framework for incorporating evolu- 
tionary processes into conservation planning using three amphib- 
ian species native to the Iberian Peninsula (Hyla molleri Bedriaga, 


1889; Pelobates cultripes (Cuvier, 1829) and Rana iberica, Boulenger, 


potential local adaptations by existing protected areas. Moreover, it can identify 
priority areas to improve conservation of evolutionary processes. Since neglecting 
evolutionary processes can impair conservation plans, we recommend using evo- 


lutionary data to inform decision-making where possible. 


adaptation, climate change, evolution, genetic diversity, prioritization, protected areas, 


resilience, single-nucleotide polymorphism 


1879). By mobilizing environmental, genetic and occurrence data, we 
mapped a comprehensive range of evolutionary attributes to delin- 
eate places containing (neutral attributes) individuals with moderate 
to high heterozygosity, different neutral genetic clusters, (adaptive 
attributes) different adaptive genetic clusters and climatic refugia. 
We then overlaid these maps with boundaries of existing protected 
areas to quantify shortfalls and identified priority areas to address 
them. We also compared these priority areas with those from 
conventional prioritization approaches that neglect evolutionary 
processes. Although challenges remain in implementing effective 
conservation policies, recent work—within our study area—shows 
how engaging with local stakeholders and governmental agencies 
can make systematic conservation planning exercises more poli- 
cy-relevant (Pinto et al., 2019). 


2 | MATERIALS AND METHODS 
2.1 | Study area and species 


Our study area encompassed the Iberian Peninsula (Figure 1). We 
examined three amphibian species with different ecological niches 
and natural histories that are endemic—or very nearly—to this area: 
the Iberian tree frog (H. molleri), western spadefoot toad (P. cultripes) 
and Iberian frog (R. iberica, Loureiro, Ferrand de Almeida, Carretero, 
& Paulo, 2008; Pleguezuelos, Marquez, & Lizana, 2002). Hyla molleri 
inhabits forests, grasslands and shrublands throughout south-west- 
ern France and most of the Iberian Peninsula, excepting eastern and 
southern Spain, and reproduces in water bodies with abundant veg- 
etation. Pelobates cultripes occurs across the Mediterranean bioregion 
of the Iberian Peninsula, south-eastern and south-western France, 
with scattered populations in the Atlantic bioregion of north-west- 
ern Iberia. It prefers soft sandy soils for burrowing, and reproduces 
in stagnant water bodies. Rana iberica occurs across northern and 
north-western Iberia, and the Sistema Central mountains in central 
Iberia. It inhabits streams and rivers with riparian vegetation, and 
small shallow water bodies. All of these species have spatially struc- 
tured patterns of neutral genetic variation (Gutiérrez-Rodriguez, 
Barbosa, & Martinez-Solano, 2017; Sanchez-Montes, Recuero, 
Barbosa, & Martinez-Solano, 2019; Teixeira, Goncalves, Ferrand, 
Garcia-Paris, & Recuero, 2018). Furthermore, R. iberica and P. cultripes 


are Near Threatened on the Red List by the International Union for 
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FIGURE 1 Panels (a-c) show maps delineating the (dark grey) 
spatial distribution of each study species within the study area and 
(points) sampling sites 


Conservation of Nature, requiring site management and protection 
(IUCN, 2019). 


2.2 | Data 


Tail and toe biopsy samples were collected from tadpoles and adults 


(respectively) from multiple sites within climatically distinct zones 


inside each species' distribution across the Iberian Peninsula and 
south-eastern France (number of samples for H. molleri: 92, P. cul- 
tripes: 184 and R. iberica: 194; Figure 1). Additional samples were ob- 
tained from museum collections (Gutiérrez-Rodriguez et al., 2017). 
ExtractMe Genomic DNA 96-Well (DNA GDANSK) and QlAamp 
DNA Micro (QIAGEN GmbtH) kits were used to extract genomic 
DNA. Depending on the amount of starting material, DNA samples 
were then concentrated with vacuum centrifugation and standard- 
ized to 300-500 ng per sample. These samples were genotyped 
by Diversity Arrays Technology Pty. Ltd, following the standard 
DArTSeq protocol (Sansaloni et al., 2011), using single-nucleotide 
polymorphisms (SNPs). The genomic datasets were then filtered 
to exclude unreliable data (see Appendix S1, Table S1, retaining a 
subset of samples and loci for H. molleri: 91 and 11,390, P. cultripes: 
178 and 25,681, and R. iberica: 169 and 24,458 respectively; 1-8 
samples per locality per species). We excluded three P. cultripes 
samples (not reported in previous statistics) because they were 
determined to have incorrect sampling locality information after 
cross-referencing preliminary results with previous work (Gutiérrez- 
Rodriguez et al., 2017). The majority of processing was completed 
using the R statistical computing environment (version 3.5.3; R Core 
Team, 2019). Genetic data were processed using PGDSpider (version 
2.1.1.5) andthe dartR R package (Gruber & Georges, 2018; Lischer 
& Excoffier, 2012). 

We compiled datasets to map evolutionary attributes for each 
species. To standardize spatial analyses and define the planning 
units for assessments and prioritizations, we created a grid with 
10 x 10 km cells covering the study area (UTM Zone 30N). Species 
distribution data were obtained from national atlases (10 x 10 km 
resolution; Figure 1; Loureiro et al., 2008; Pleguezuelos et al., 2002). 
Seven bioclimatic variables were obtained to characterize contempo- 
rary climate regimes (1950-1990; BIO1: annual mean temperature, 
BIO2: mean diurnal range, BIO3: isothermality, BIO8: mean tempera- 
ture of wettest quarter, BIO9: mean temperature of driest quarter, 
BIO13: precipitation of wettest month and BlO14: precipitation of 
driest month; 2.5’ resolution; Hijmans, Cameron, Parra, Jones, & 
Jarvis, 2005). These variables were selected from the 19 available 
bioclimatic variables to avoid multicollinearity issues (following 
Naimi, Hamm, Groen, Skidmore, & Toxopeus, 2014). Projections of 
the seven bioclimatic variables from 19 general circulation models 
were also obtained to characterize plausible future climatic regimes 
(2070), assuming representative concentration pathway 4.5 (2.5' res- 
olution; Hijmans et al., 2005). Additionally, soil bedrock data were 
obtained to characterize abiotic conditions influencing vegetation 
communities (1 km? resolution; Panagos, 2006; Van Liedekerke, 
Jones, & Panagos, 2006). Since 3.94% of the planning units (grid 
cells) were missing bedrock classification data, missing data were 
interpolated using modal values of adjacent locations. The datasets 
were then reprojected and aligned to the spatial grid (Figures S1- 
S20). Spatial processing was performed using the sf and raster R 
packages (Hijmans, 2018; Pebesma, 2018). 

We also compiled datasets to assess existing protected areas and 


generate prioritizations. Data delineating the boundaries of protected 
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areas in the study area were also obtained (IUCN & UNEP-WCMC, 
2018) and cleaned following standard practices using the wdpar R 
package (Appendix S2; Hanson, 2019). Additionally, 1,574 Sites of 
Community Importance (Habitat Directive), 709 Special Protection 
Areas (Birds Directive) and 285 Natura 2000 sites were excluded 
because such places are not strictly managed as wildlife refugia (e.g. 
hunting is permitted; DG XI.D.2, 1998). The most recent version of 
the human footprint index, describing anthropogenic pressure during 
2009, was obtained to characterize opportunity costs for establishing 
new reserves (1 km? resolution; Venter et al., 2016, 2018). Since 1.23% 
of the study area were missing human footprint data, an inverse dis- 
tance weighting interpolation procedure was used to predict missing 
values. These datasets were then reprojected and aligned to the spatial 
grid, and grid cells with over 50% coverage by protected areas were 


subsequently treated as protected (Figures S21 and $22). 


2.3 | Statistical analysis 


Many species comprise multiple distinct evolutionary units (Funk 
et al., 2012). To account for such units, we grouped the individuals 
sampled for each species into broad-scale genetic clusters (here- 
after, genetic lineages) using Structure (version 2.3.4; Pritchard, 
Stephens, & Donnelly, 2000), clumpp (version 1.1.2; Jakobsson 
& Rosenberg, 2007) and Evanno's AK method (Evanno, Regnaut, 
& Goudet, 2005) implemented in the pophelper R package 
(Appendix S3 and Figures $23-S24; Francis, 2017). We also veri- 
fied these analyses using the tess3r R package (Figure S25; 
Caye, Jay, Michel, & Francois, 2018). Next, we used phylogeo- 
graphic interpolation to predict the spatial distribution of each 
genetic lineage—per the Structure analyses—for each species 
using the phylin R package (Tarroso, Carvalho, & Velo-Anton, 
2019). For each species, semi-variograms were estimated using 
exponential, Gaussian, linear and spherical models fitted with a 
maximum of 100,000 iterations and evaluated using coefficient 
of determination statistics (Table $2). The best semi-variogram fit 
for each species was then used with ordinary kriging to predict 
the spatial distribution of each genetic lineage for each species 
(Figures S26-S27). 

We identified loci potentially under selection using outlier loci 
detection and environmental association analyses (Appendix $4). 
This step was completed after grouping the samples into genetic 
lineages because one of these analyses requires such information. 
To prevent individuals of uncertain ancestry from biasing these 
analyses, they were excluded (i.e. individuals with a maximum 
membership probability below 75%). The following environmental 
association analyses were performed using the contemporary cli- 
mate data: latent factor mixed models (using the 1fmm R package; 
Caye, Jumentier, Lepeule, & Francois, 2019), and redundancy anal- 
yses (Forester, Lasky, Wagner, & Urban, 2018). The following out- 
lier detection analyses were also performed: multinomial-Dirichlet 
models (via BayeScan version 2.1, and the qvalue R package; Foll 
& Gaggiotti, 2008; Storey, Bass, Dabney, & Robinson, 2015), Moran 


spectral outlier detection analyses (using the spdep R pack- 
age; Bivand, Pebesma, Gomez-Rubio, & Pebesma, 2013; Wagner, 
Chavez-Pesqueira, & Forester, 2017) and principal components 
analyses (using the pcadapt, nFactor R,and qvalue R packages; 
Luu, Bazin, & Blum, 2017; Raiche, 2010). After running these analy- 
ses (Figures S28-S32), loci flagged in at least two of the outlier loci 
analyses or two of the environmental association analyses were 
subsequently treated as putatively adaptive and the remainder as 
neutral (number of adaptive loci for H. molleri: 18, P. cultripes: 151 
and R. iberica: 215; Table S3). 

We fitted environmental niche models for each species using 
MaxEnt (Hirzel & Le Lay, 2008; Phillips, Dudík, & Schapire, 2017). 
After fitting the models using the atlas and contemporary environ- 
mental data (using the ENMeval R package; Muscarella et al., 2014), 
they were used to map contemporary habitat suitability for each 
species (Appendix S5, Figure S33; Table $4). They were also used 
to produce 19 maps of predicted future habitat suitability for each 
species according to the future bioclimatic projections and contem- 
porary bedrock data (Figures $34-S36). 


2.4 | Evolutionary attributes 


We mapped four evolutionary attributes for each species. 
Specifically, two evolutionary attributes related to neutral processes 
(individual heterozygosity and neutral genetic clusters) and two evo- 
lutionary attributes related to adaptive processes (adaptive genetic 


clusters and climatic refugia). See Table 1 for rationale. 


2.4.1 | Individual heterozygosity 


For each species, we (a) calculated heterozygosity scores for each 
individual using neutral loci (using the inbreedR R package; 
Coltman, Pilkington, Smith, & Pemberton, 1999; Stoffel et al., 2016), 
(b) computed the average heterozygosity of each sampling locality 
to facilitate spatial interpolation, (c) interpolated the heterozygosity 
scores for each planning unit inside the species' distribution using 
thin plate splines (Figure S37; Table S5) and (d) for each genetic line- 
age, clamped values for planning units with heterozygosity scores 
below the median value to zero to exclude planning units with low 


heterozygosity (Figure $38). 


2.4.2 | Neutral genetic clusters 


For each genetic lineage within each species, we (a) used hier- 
archical genetic cluster analyses to classify individuals into fine- 
scale genetic clusters (hereafter, genetic clusters) using neutral 
loci (using Gaussian mixture model-based cluster analyses im- 
plemented in the mclust R package; Scrucca, Fop, Murphy, & 
Raftery, 2016; following Van den Wyngaert, Most, Freimann, 
Ibelings, & Spaak, 2015), (b) estimated semi-variograms for 


HANSON ET AL. 


Journal of Applied Ecology 


5 





TABLE 1 Description of evolutionary attributes 


Attribute Rationale 


Neutral evolutionary processes 


Individual heterozygosity By conserving planning units that contain individuals with 
moderate to high heterozygosity, protected area systems 
may secure individuals with high genetic fitness that 
belong to demographically stable populations (Grueber 


et al., 2008) 


Neutral genetic clusters By conserving a set of planning units that contain 
individuals with different combinations of alleles, 
protected area systems can limit erosion of overall 


genetic diversity (Carvalho et al., 2011; Moritz, 2002) 
Adaptive evolutionary processes 


Adaptive genetic clusters By conserving a set of planning units that contain 
individuals with different combinations of alleles at loci 
that are (likely) under selection, protected area systems 
can prevent loss of locally adaptive genetic variants 


(Sgro et al., 2011) 


Climate refugia By conserving planning units that are expected to 
experience relatively slow declines in environmental 
suitability, protected area systems may secure places 


that provide more time for species to adapt to changing 


Targets 


Represent 15% of the planning units with moderate 
to high individual heterozygosity scores, weighted 
by the scores, associated with each genetic lineage 
in each species (Figure S38) 


Represent 15% of the planning units associated 
with each genetic cluster in each species 
(Figures S39-S41) 


Represent 15% of the planning units associated 
with each genetic cluster, weighted by 
the non-neutrality scores, in each species 
(Figures S42-S44) 


Represent 15% of the planning units with climatic 
refugia, weighted by long-term suitability 
scores, for each genetic lineage in each species 
(Figure S46) 


conditions (Game, Wallis, & Jamieson, 2011) 


phylogeographic interpolation following the previously described 
methodology (Table S6), (c) predicted the spatial distribution of 
each genetic cluster within the lineage (using the phylin R pack- 
age) and (d) for each species, associated each planning unit with 
the most likely genetic cluster to occupy it (using cluster assign- 
ments for sampled planning units and interpolated assignments 
for those missing samples; Figures S39-S41). We used Gaussian 
mixture model-based cluster analyses since they can be applied to 


both neutral and adaptive loci. 


2.4.3 | Adaptive genetic clusters 


For each species, we employed the same methodology described 
above—except using putatively adaptive loci-to map genetic 
clusters within each genetic lineage (Table S6). To account for 
spatial variation in the strength of adaptive processes: for each 
species, we (a) performed metric multi-dimensional scaling analy- 
ses with Manhattan distances (Table S7) to produce separate 
ordinations for the putatively adaptive and neutral loci (using 
the vegan R package; Oksanen et al., 2018), (b) subjected the 
ordinations to a Procrustes analysis, (c) computed the absolute 
values of the residuals from the Procrustes analysis to create 
non-neutrality scores that assign greater values to individuals 
that contain combinations of alleles at putatively adaptive loci 
that deviate more strongly from patterns among neutral loci, (d) 
spatially interpolated the scores across each species’ distribution 
using thin plate splines (Table S8) and (e) used the non-neutrality 
scores to weight planning units within each adaptive genetic clus- 
ter (Figures S42-S44). 


2.4.4 | Climatic refugia 


For each species, we combined maps of contemporary and potential 
future habitat suitability (Figures S33-S36) to produce a single map 
denoting long-term habitat suitability (Figure S45). Specifically, we 
used harmonic means because they penalize localities with high vari- 
ation in suitability scores. Thus, places containing higher long-term 
habitat suitability scores are more likely to serve as future climatic 
refugia. For each genetic lineage, we clamped planning units with 
scores below the median value to zero to exclude planning units with 


less suitable conditions (Figure $46). 


2.5 | Protected area assessments and prioritizations 


We assessed how well the evolutionary attributes are covered by 
existing protected areas (termed representation) and compared their 
percent coverage to target thresholds to reveal shortfalls. Ideally, 
targets would be derived using mathematical models related to 
species' persistence. However, no such models currently exist and 
so, following previous work (e.g. Carvalho et al., 2011; Thomassen 
et al., 2011), we set fixed targets to represent 15% of each evolution- 
ary attribute for each species (see Table 1). 

We generated a prioritization using the evolutionary attributes 
(using the prioritizr R package; Hanson, Schuster, et al., 2019). 
To assess the relative importance of priority areas, we derived irre- 
placeability scores using a modified version of the replacement cost 
method (Appendix S6; Cabeza & Moilanen, 2006). To gain insight into 
their roles, we identified which evolutionary attribute was helped the 


most by each priority area (measured using percentage-wise increase 
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FIGURE 2 Maps show (a) priority areas for representing the species' evolutionary attributes, (b) their irreplaceability scores, (c) the 
primary evolutionary attribute for which they increase representation and (d) the prioritization generated following conventional approaches 
that neglect evolutionary processes. Labels in (c) correspond to different regions of priority areas, and the line demarcates regions 1 and 2 


in representation). We visually grouped the priority areas primarily as- 
sociated with the same evolutionary attributes into regions to facilitate 
interpretation (shown in Figure 2c). To understand how well conven- 
tional approaches represent evolutionary processes, we generated a 
prioritization using only the species' distribution data—omitting the 
evolutionary attributes—to secure 15% of each species' distribution 
within the study area. All prioritizations were generated using the min- 
imum set formulation of the reserve selection problem with planning 
unit costs set as the human footprint index and solved to optimality 
(using Gurobi version 8.1.0; Gurobi Optimization, LLC, 2018). 


3 | RESULTS 


Existing protected areas in the Iberian peninsula are failing to ad- 
equately represent evolutionary attributes for H. molleri, P. cultripes 


and R. iberica (Table 2; Table S8). Each of these species comprise two 


distinct genetic lineages. Protected areas are not adequately repre- 
senting places predicted to contain individuals with moderate to high 
heterozygosity—none of the genetic lineages for any of the three 
species have adequate representation of such places. When examin- 
ing fine-scale patterns of intraspecific genetic variation within each 
of the species' genetic lineages, less than one-third of the species' 
genetic clusters classified using neutral loci are adequately repre- 
sented (H. molleri: 16.67%, P. cultripes: 12.5% and R. iberica: 25%). 
Additionally, none of the genetic clusters identified using putatively 
adaptive loci are adequately represented for H. molleri and P. cul- 
tripes, and only half of such genetic clusters are adequately repre- 
sented for R. iberica. Furthermore, none of the genetic lineages for 
H. molleri and P. cultripes and only one of the two genetic lineages 
comprising R. iberica have adequately represented climatic refugia. 
Priority areas were identified for improving the representation 
of the species' evolutionary attributes by protected areas (Figure 2a; 


Table 2; Table S9). If added to the existing protected area system, 
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TABLE 2 Performance of the existing protected area system, prioritization generated using evolutionary attributes and conventional 
prioritization that neglected evolutionary processes. For a given species, each evolutionary attribute is associated with multiple genetic 

lineages (identified using the Structure software) or multiple genetic clusters (identified using the mclust R package). Data show the 
number and percentage of clusters or lineages that are adequately represented by a given conservation plan for a given evolutionary 


attribute. Bold font indicates adequate representation 


Species Process Attribute 
Hyla molleri Neutral Individual heterozygosity 
Neutral genetic clusters 
Adaptive Adaptive genetic clusters 
Climatic refugia 
Pelobates Neutral Individual heterozygosity 
cultripes Neutral genetic clusters 
Adaptive Adaptive genetic clusters 
Climatic refugia 
Rana iberica Neutral Individual heterozygosity 
Neutral genetic clusters 
Adaptive Adaptive genetic clusters 


Climatic refugia 


these places would increase the size of the protected area estate 
from 8.65% to 12.1% of the Iberian Peninsula. Most priority areas 
have relatively low irreplaceability scores (<3; Figure 2b) and so they 
can be substituted by places with low human pressure. Many of the 
priority areas are sited at the intersection of the northern genetic 
lineages of H. molleri and R. iberica (region 1, Figure 2c). These spe- 
cific locations—in addition to helping to address shortfalls for other 
evolutionary attributes—markedly improve the capacity of existing 
protected areas to represent climatic refugia and different adap- 
tive and neutral genetic clusters. Moving south, more priority areas 
are sited at the intersection of the southern lineages of these two 
species (region 2, Figure 2c). They mainly help increase the repre- 
sentation of adaptive and neutral genetic clusters, and individuals 
with moderate to high heterozygosity. Further east, more priority 
areas are sited in places containing the eastern genetic lineage for P. 
cultripes (region 3, Figure 2c). They primarily help increase the repre- 
sentation of climatic refugia and adaptive and neutral genetic clus- 
ters (similar to region 1). These results suggest that different priority 
areas may be important for different evolutionary processes. 

The prioritization generated using conventional approaches that 
neglect evolutionary processes had poor coverage of the species’ evo- 
lutionary attributes (Figure 2d; Table 2; Table $10). It only adequately 
represented places containing individuals with moderate to high het- 
erozygosity for half of the genetic lineages comprising H. molleri and 
P. cultripes and none of the genetic lineages comprising R. iberica. It 
was also only able to adequately represent half of the neutral genetic 
clusters associated with each of the three species. Additionally, it only 
adequately represented one-third of the adaptive genetic clusters for 
H. molleri and P. cultripes (both 33.33%), and provided no improvement 
for R. iberica. Similarly, it increased coverage to adequately represent 


climatic refugia for half of the genetic lineages comprising H. molleri 


Number of Existing 

clusters/ protected Evolutionary Conventional 
lineages areas prioritization prioritization 
2 0 (0%) 2 (100%) 1 (50%) 

6 1 (16.67%) 6 (100%) 3 (50%) 

9 0 (0%) 9 (100%) 3 (33.33%) 

2 0 (0%) 2 (100%) 1 (50%) 

2 0 (0%) 2 (100%) 1 (50%) 

8 1 (12.5%) 8 (100%) 4 (50%) 

9 0 (0%) 9 (100%) 3 (33.33%) 

2 0 (0%) 2 (100%) 1 (50%) 

2 0 (0%) 2 (100%) 0 (0%) 

8 2 (25%) 8 (100%) 4 (50%) 

6 3 (50%) 6 (100%) 3 (50%) 

2 1 (50%) 2 (100%) 1 (50%) 


and P. cultripes, and provided no improvement for R. iberica. Despite 
performing much worse than the prioritization generated using evolu- 
tionary attributes, this prioritization encompassed a similar percentage 
of the study area (11.16%) and cost (i.e. sum human footprint index 
values) only 8.88% less. These findings demonstrate that prioritiza- 
tions generated without explicitly considering evolutionary processes 


can fail to represent them. 


4 | DISCUSSION 
We developed a new framework for incorporating adaptive and 
neutral evolutionary processes into systematic conservation 
planning and applied it to three amphibian species (H. molleri, P. 
cultripes and R. iberica). Although existing protected areas in our 
study area are helping to conserve biodiversity (Araujo, Lobo, 
& Moreno, 2007), we found that additional protected areas are 
needed to help maintain adaptive and neutral evolutionary pro- 
cesses for these species. We reveal priority areas that would 
achieve adequate representation of such places for only a minor 
increase in the protected area system. We also show that ne- 
glecting evolutionary processes can produce prioritizations that 
fail to substantially improve the representation of evolutionary 
processes—further highlighting the importance of accounting for 
them (Carvalho et al., 2011; Hanson, Fuller, et al., 2019). Our work 
helps pave the way for the widespread integration of evolution- 
ary data into conservation decisions (Beger et al., 2014; Carvalho 
et al., 2017; Cowling et al., 2003). 

Our framework enables conservation planners to explicitly and 
comprehensively account for adaptive and neutral evolutionary pro- 


cesses. We built this framework on the strengths and weaknesses of 
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previous frameworks (e.g. Beger et al., 2014; Diniz-Filho et al., 2012; 
Hanson et al., 2018). It accounts for the evolutionary processes op- 
erating within each species separately (Hanson et al., 2018). Instead 
of relying on surrogates (Carvalho et al., 2011; Ponce-Reyes, Clegg, 
Carvalho, McDonald-Madden, & Possingham, 2014), it uses metrics 
derived from putatively adaptive and neutral loci (Hanson et al., 2017). 
Additionally, it uses genetic clusters—instead of allele frequencies 
(Diniz-Filho et al., 2012)—to characterize fine-scale patterns of genetic 
diversity so that spatial interpolation routines can be used to make 
predictions for planning units missing data (Carvalho et al., 2017). It 
also emphasizes places containing individuals with greater genetic di- 
versity (Beger et al., 2014; Nielsen et al., 2017; Paz-Vinas et al., 2018; 
Thomassen et al., 2011). It accounts for abiotic factors that influ- 
ence evolutionary processes (i.e. climatic refugia; Game et al., 2011). 
Furthermore, by leveraging the minimum set formulation of the re- 
serve selection problem (Beger et al., 2014; Nielsen et al., 2017), opti- 
mal solutions can be obtained quickly for planning exercises involving 
thousands of planning units (unlike Hanson et al., 2018). 

We found that existing protected areas in the Iberian Peninsula 
are not adequately representing evolutionary attributes for the 
three amphibian species. This result is not surprising given that 
the existing protected area system was not designed specifically 
for these species—much less their evolutionary processes. Many 
protected areas in Europe are in places that historically served as 
royal hunting grounds, gardens or forests (European Environment 
Agency, 2012). Additionally, they may also have been established 
opportunistically in places that provide limited benefits (Joppa & 
Pfaff, 2009). We also identified priority areas to improve the rep- 
resentation of the species' evolutionary attributes. Since we ex- 
amined far fewer species than previous studies—and differences 
in the spatial patterns of intraspecific genetic diversity among 
species can alter priorities (Paz-Vinas et al., 2018)—few of the pri- 
ority areas identified in this study coincide precisely with those in 
previous studies (Araújo et al., 2007; Carvalho et al., 2011, 2017). 
Despite such fine-scale differences, our study reinforces previ- 
ous work that has identified central Portugal, and north-western 
(Galicia) and eastern (Aragón and Comunidad Valenciana) Spain 
as important conservation areas in the Iberian Peninsula (Araújo 
et al., 2007; Carvalho et al., 2017). 

Our results show that conventional prioritization approaches that 
do not explicitly incorporate evolutionary attributes may fail to main- 
tain evolutionary processes. These results contribute to the growing 
body of evidence indicating that neglecting evolutionary processes 
can undermine conservation (Carvalho et al., 2011; Hanson, Fuller, 
et al., 2019; Nielsen et al., 2017; Thomassen et al., 2011). For in- 
stance, within the Iberian Peninsula, prioritizations generated using 
only species’ distribution data secure much less intraspecific neu- 
tral genetic diversity and environmental heterogeneity—a putative 
surrogate for adaptive genetic diversity (Hanson et al., 2017)—than 
prioritizations generated using these evolutionary attributes directly 
(Carvalho et al., 2011, 2017). Outside of the Iberian Peninsula, pri- 
oritizations generated using such conventional approaches secure 


much less intraspecific genetic diversity and morphological variation 


(Nielsen et al., 2017; Thomassen et al., 2011; but see Hermoso 
et al., 2016). Furthermore, priority areas identified at the global-scale 
using such conventional approaches secure much less of the diver- 
sity in environmental conditions found across species’ distributions 
(Hanson, Rhodes, et al., 2020)—meaning they could omit important 
places for climate change adaptation (Sgro et al., 2011). 

Our study has limitations. First, while real-world planning ex- 
ercises should consider far more than three species (Margules & 
Pressey, 2000), our main aim was to introduce the framework—not 
to identify priority areas for conserving biodiversity across the 
Iberian Peninsula. Second, one remaining challenge is identifying 
effective targets for evolutionary attributes that will substantially 
improve species' long-term persistence. Although we adopted fixed 
targets following previous studies (e.g. Carvalho et al., 2011), differ- 
ent targets will likely be needed for different evolutionary attributes 
(Pressey, Cowling, & Rouget, 2003). Third, although samples were 
collected to maximize our capacity to map adaptive genetic varia- 
tion, we were unable to assess population-level genetic diversity (e.g. 
allelic richness). Fourth, our five analyses for detecting putatively 
adaptive loci did not yield highly consistent results (further explored 
in A. Marques, J.O. Hanson, M. Camacho-Sanchez, Í. Martínez- 
Solano, C. Moritz; P. Tarroso, G. Velo-Antón; A. Verissímo, & S.B. 
Carvalho, unpubl. data), and so further experimentation is required. 
Additionally, we detected far fewer putatively adaptive loci for H. 
molleri, and so our results may not reflect its full range of adaptive 
genetic variation. Fifth, although we followed standard practices for 
modelling future species' distributions, uncertainty in future climate 
projections limits our ability to correctly identify climatic refugia 
(Game et al., 2011). Sixth, establishing protected areas alone is often 
insufficient, and populations may also require direct management in- 
terventions (e.g. genetic rescue; Flanagan, Forester, Latch, Aitken, & 
Hoban, 2018). Finally, future research could examine gene flow (e.g. 
using Daigle et al., 2020; Hanson, Fuller, et al., 2019). 

Conservation plans need to maintain adaptive and neutral evolu- 
tionary processes to maximize the persistence of biodiversity. Since 
challenges remain in explicitly using these processes to inform de- 
cision-making, we developed a framework to identify priority areas 
using ecological and evolutionary data. Further research is needed 
to inform setting representation targets for evolutionary processes. 
In the meantime, we advise conservation scientists and practitioners 
to carefully consider using evolutionary processes to guide reserve 


selection. 
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